!>\brief
!!
!! Harmonic map equation.
!!
!! \n
!!
!! We consider the problem
!!
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!  \frac{d}{dt}\int_K \underline{d}\cdot\bar{\underline{d}}\,dx +
!!  \int_K \nabla\underline{d}:\nabla\bar{\underline{d}}\,dx -
!!  \int_{\partial K}
!!  \left\{\nabla\underline{d}\right\}:\llbracket\bar{\underline{d}}\rrbracket\,d\sigma
!!  -\kappa \int_{\partial K}
!!  \left\{\nabla\bar{\underline{d}}\right\}:\llbracket\underline{d}\rrbracket\,
!!  d \sigma \\[4mm]
!!  \displaystyle
!!  +\int_{\partial K}
!!  \alpha\llbracket\underline{d}\rrbracket:\llbracket\bar{\underline{d}}\rrbracket\,d\sigma
!!  +\int_K q\underline{d}\cdot\bar{\underline{d}}\,dx & = & 0
!!  \\[4mm]
!!  \displaystyle
!!  \int_K d^2\bar{q}\,dx & = & \displaystyle \int_K\bar{q}\,dx
!!  \end{array}
!! \f}
!!
!! Due to the presence of the Lagrange multiplier, we assume here that
!! an implicit solver is used.
!!
!! Enforcing the boundary conditions in weak form, the above system
!! can be reduced to the DAE
!!
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!  M \dot{\tt d} + A{\tt d} + B^T_{\tt d}{\tt q} & = &
!!  \displaystyle {\tt f} \\[2mm]
!!  \displaystyle B_{\tt d}{\tt d} & = & {\tt u}
!!  \end{array}
!! \f}
!! To recast this in the form required by the BDF solvers, we write
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!  \dot{\tt d} & = & \displaystyle - M^{-1}A{\tt d} - M^{-1}B^T_{\tt
!!  d}{\tt q} + M^{-1}{\tt f} \\
!!  \displaystyle B_{\tt d}{\tt d} - {\tt u}& = & 0
!!  \end{array}
!! \f}
!! so that we have
!! \f{displaymath}{
!!  f({\tt d},{\tt q}) = - M^{-1}A{\tt d} - M^{-1}B^T_{\tt
!!  d}{\tt q} + M^{-1}{\tt f}, \qquad
!!  g({\tt d}) = B_{\tt d}{\tt d} - {\tt u}.
!! \f}
!! Hence, at each time step we need to solve
!! \f{displaymath}{
!!  \begin{array}{rcl}
!!  \displaystyle
!!  M{\tt d} + \sigma A{\tt d} + \sigma B^T_{\tt d}{\tt q} & = &
!!   {\tt b}_{\tt d} + \sigma{\tt f} \\[2mm]
!!  \displaystyle B_{\tt d}{\tt d} & = & {\tt u}
!!  \end{array}
!! \f}
!! or, in matrix form,
!! \f{displaymath}{
!!  \left[
!!  \begin{array}{cc}
!!  M + \sigma A & \sigma B^T_{\tt d} \\
!!  \sigma B_{\tt d} & 0
!!  \end{array}
!!  \right]
!!  \left[
!!  \begin{array}{c}
!!   {\tt d} \\
!!   {\tt q}
!!  \end{array}
!!  \right] =
!!  \left[
!!  \begin{array}{c}
!!   {\tt b}_{\tt d} + \sigma{\tt f} \\
!!   \sigma{\tt u} 
!!  \end{array}
!!  \right].
!! \f}
!<----------------------------------------------------------------------
module mod_hmap_ode

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp
 
 use mod_linal, only: &
   mod_linal_initialized, &
   linsys

 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_hmap_state, only: &
   mod_hmap_state_initialized, &
   t_hmap_state, source
 
 use mod_time_integrators, only: &
   mod_time_integrators_initialized, &
   c_ode

 use mod_sparse, only: &
   mod_sparse_initialized!, &
 !  ! sparse types
 !  t_col,       &
 !  t_tri,       &
 !  ! construction of new objects
 !  new_col,     &
 !  new_tri,     &
 !  ! convertions
 !  col2tri,     &
 !  tri2col,     &
 !  tri2col_skeleton, &
 !  tri2col_skeleton_part, &
 !  ! overloaded operators
 !  operator(+), &
 !  operator(*), &
 !  sum,         &
 !  transpose,   &
 !  matmul,      &
 !  ! error codes
 !  wrong_n,     &
 !  wrong_m,     &
 !  wrong_nz,    &
 !  wrong_dim,   &
 !  ! other functions
 !  nnz_col,     &
 !  nz_col,      &
 !  nz_col_i,    &
 !  get,         &
 !  set,         &
 !  diag,        &
 !  spdiag,      &
 !  ! deallocate
 !  clear

 !use mod_octave_io_sparse, only: &
 !  mod_octave_io_sparse_initialized, &
 !  write_octave

 !use mod_linsolver, only: &
 !  mod_linsolver_initialized, &
 !  c_linpb, c_itpb, c_mumpspb, c_umfpackpb, &
 !  gmres

 use mod_newton, only: &
   mod_newton_initialized, &
   t_nwt_params, c_nwt, &
   newton

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_hmap_ode_constructor, &
   mod_hmap_ode_destructor,  &
   mod_hmap_ode_initialized, &
   t_hmap_ode

 private

!-----------------------------------------------------------------------

! Module types and parameters

 !> Harmonic map ODE
 type, extends(c_ode) :: t_hmap_ode
 contains
  procedure, nopass :: rhs   => hmap_rhs
  procedure, nopass :: solve => hmap_solve
  procedure, nopass :: scal  => hmap_scal
 end type t_hmap_ode

 !> Newton method used at each iteration
 type, extends(c_nwt) :: t_hmap_nwt
  real(wp) :: res(3)
  logical :: residual_set = .false.
 contains
  procedure, pass(z) :: pres => nwt_pres
  procedure, pass(x) :: nres => nwt_nres
  procedure, pass(x) :: norm => nwt_norm
  procedure, pass(y) :: defy => nwt_defy
  procedure, pass(x) :: source => nwt_source
  !> Update the field \c res to be consistent with \c x
  procedure, pass(x) :: update_res
 end type t_hmap_nwt
 
! Module variables

 !> At each call to <tt>t_hmap_ode\%solve</tt>, one must solve a
 !! nonlinear problem. This type collects the parameters wich define
 !! such problem, so that they are available to the type bound
 !! procedures included in \c t_hmap_nwt.
 type t_newton_data
  real(wp) :: sigma
  real(wp) :: rhs(3)
 end type t_newton_data
 type(t_newton_data) :: newton_data

 ! Linear system
 type(t_col) :: mmmt !< transposed matrix of the linear system
 integer, allocatable :: tom(:) !< indexes to mmmt%x

 logical, protected ::               &
   mod_hmap_ode_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_hmap_ode'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_hmap_ode_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
  (mod_state_vars_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
  (mod_hmap_state_initialized.eqv..false.) .or. &
(mod_time_integrators_initialized.eqv..false.) .or. &
      (mod_sparse_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_hmap_ode_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   !--------------------------------------------------------------------
   ! Define the matrix sparsity pattern
   d  = grid%d
   pk = base%pk
   pk2 = pk**2
   elem_do: do ie=1,grid%ne

     ! mass and gradient matrices
     do j=1,pk; do jd=1,d
       do i=1,pk; do id=1,d
         pos =  (ie-1)*(pk**2)
         iii(pos) = (i-1)*d + id
         jjj(pos) = (j-1)*d + jd
       enddo
     enddo
      


   enddo elem_do
   !--------------------------------------------------------------------


   mod_hmap_ode_initialized = .true.
 end subroutine mod_hmap_ode_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_hmap_ode_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_hmap_ode_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_hmap_ode_initialized = .false.
 end subroutine mod_hmap_ode_destructor

!-----------------------------------------------------------------------

 !> RHS not implemented: use \c hmap_solve instead.
 subroutine hmap_rhs(tnd,t,uuu)
  real(wp),     intent(in)    :: t   !< time level
  class(c_stv), intent(in)    :: uuu !< present state
  class(c_stv), intent(inout) :: tnd !< tendency

 end subroutine hmap_rhs

!-----------------------------------------------------------------------

 subroutine hmap_solve(t,x,sigma,b,xl)
  real(wp),     intent(in) :: t, sigma
  class(c_stv), intent(in) :: b, xl
  class(c_stv), intent(inout) :: x

  type(t_hmap_nwt) :: xx
  type(t_nwt_params) :: params
  logical :: lmaxiter
  integer :: niter
  real(wp), allocatable :: res(:,:)

  real(wp) :: m(3,3), sol(3)

   params%nmax = 50
   params%tolerance = 1e-6_wp

!   select type(x); type is(t_hmap_state)
!     select type(b); type is(t_hmap_state)
!       select type(xl); type is(t_hmap_state)
!
!   m(1,:) = (/ 1.0_wp ,  sigma , sigma*xl%d(1,1) /)
!   m(2,:) = (/ -sigma , 1.0_wp , sigma*xl%d(2,1) /)
!   m(3,:) = (/ xl%d(1,1) , xl%d(2,1) , 0.0_wp /)
!
!   newton_data%rhs(1:2) = b%d(:,1)
!   newton_data%rhs(3) = 1.0_wp
!
!   call linsys(m,newton_data%rhs,sol)
!
!   x%d(:,1) = sol(1:2)
!   x%q = sol(3)
!
!       end select
!     end select
!   end select

   select type(b); type is(t_hmap_state)

   ! Set the system variables defining the nonlinear problem
   newton_data%sigma = sigma
   newton_data%rhs = (/ b%d(:,1) , 1.0_wp /)
   ! Provide a meaningful initial guess
   call x%copy(xl)
   ! Solve the nonlinear problem
   call newton(x,xx, params, lmaxiter, niter, res)

   write(*,*) 'Solving for time ',t,'; iterations ',niter
   write(*,*) 'residuals: ',res(1,:)

   end select

 end subroutine hmap_solve
 
!-----------------------------------------------------------------------

 function hmap_scal(x,y) result(s)
  class(c_stv), intent(in) :: x, y
  real(wp) :: s

   select type(x); type is(t_hmap_state)
     select type(y); type is(t_hmap_state)

   s = sum(x%d*y%d) + sum(x%q*y%q)

     end select
   end select
 end function hmap_scal

!-----------------------------------------------------------------------

 subroutine nwt_pres(z,x,eta,linres,liniter)
  class(c_nwt), intent(inout) :: x
  class(t_hmap_nwt), intent(inout) :: z
  real(wp), intent(in)  :: eta    !< forcing term
  real(wp), intent(out) :: linres !< linear solver residual
  integer, intent(out) :: liniter !< linear solver iterations

  real(wp) :: sigma
  real(wp) :: m(3,3), sol(3)

   ! Recover the required parameters
   sigma = newton_data%sigma

   select type(x); type is(t_hmap_nwt)
     select type(xx=>x%x); type is(t_hmap_state)
     select type(zx=>z%x); type is(t_hmap_state)

   m(1,:) = (/ 1.0_wp ,  sigma , sigma*xx%d(1,1) /)
   m(2,:) = (/ -sigma , 1.0_wp , sigma*xx%d(2,1) /)
   m(3,:) = (/ xx%d(1,1) , xx%d(2,1) , 0.0_wp /)

   if(.not.x%residual_set) call x%update_res()

   call linsys(m,-x%res,sol)

   zx%d(:,1) = sol(1:2)
   zx%q = sol(3)

   linres = 0.0_wp
   liniter = 1

     end select
     end select
   end select
 end subroutine nwt_pres

!-----------------------------------------------------------------------

 function nwt_nres(x) result(s)
  class(t_hmap_nwt), intent(inout) :: x
  real(wp) :: s

  real(wp) :: sigma

   if(.not.x%residual_set) call x%update_res()
   s = sqrt(sum(x%res**2))

 end function nwt_nres

!-----------------------------------------------------------------------

 function nwt_norm(x) result(n)
  class(t_hmap_nwt), intent(in) :: x
  real(wp) :: n

   select type(xx=>x%x); type is(t_hmap_state)

   n = sqrt( sum(xx%d**2) + sum(xx%q**2) )

   end select
 end function nwt_norm

!-----------------------------------------------------------------------

 subroutine nwt_defy(y,x)
  class(c_stv),      intent(in)    :: x
  class(t_hmap_nwt), intent(inout) :: y

   call y%x%copy(x)
   y%residual_set = .false.
 end subroutine nwt_defy

!-----------------------------------------------------------------------

 subroutine nwt_source(y,x)
  class(t_hmap_nwt), intent(in)          :: x
  class(c_nwt), allocatable, intent(out) :: y

   allocate(t_hmap_nwt::y)

   call x%x%source(y%x)

   select type(y); type is(t_hmap_nwt)
   y%residual_set = .false.
   end select

 end subroutine nwt_source

!-----------------------------------------------------------------------
 
 subroutine update_res(x)
  class(t_hmap_nwt), intent(inout) :: x

  real(wp) :: sigma
  character(len=*), parameter :: &
    this_sub_name = 'update_res'
 
   ! Recover the required parameters
   sigma = newton_data%sigma

   select type(xx=>x%x); type is(t_hmap_state)

   x%res(1) = xx%d(1,1) + sigma*xx%d(2,1) + sigma*xx%d(1,1)*xx%q(1) &
           - newton_data%rhs(1)
   x%res(2) = xx%d(2,1) - sigma*xx%d(1,1) + sigma*xx%d(2,1)*xx%q(1) &
           - newton_data%rhs(2)
   x%res(3) = xx%d(1,1)**2+xx%d(2,1)**2 - 1.0_wp

   end select

   x%residual_set = .true.

 end subroutine update_res
 
!-----------------------------------------------------------------------

end module mod_hmap_ode

